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We study synchronization phenomenon in a self-correcting population of noisy phase oscillators 
with randomly distributed natural frequencies. In our model each oscillator stochastically switches 
its phase to the ensemble-averaged value ip at a typical rate which is proportional to the degree of 
phase coherence r. The system exhibits a continuous phase transition to collective synchronization 
similar to classical Kuramoto model. Based on the self-consistent arguments and on the linear 
stability analysis of an incoherent state we derive analytically the threshold value k c of coupling 
constant corresponding to the onset of a partially synchronized state. Just above the transition 
point the linear scaling law r oc k — k c is found. We also show that nonlinear relation between 
rate of phase correction and order parameter leads to non-trivial transition between incoherence 
and synchrony. To illustrate our analytical results, numerical simulations have been performed for a 
large population of phase oscillators with proposed type of coupling. The results of this work could 
become useful in designing distributed networked systems capable of self-synchronization. 


I. INTRODUCTION 


The spontaneous synchronization of mutually coupled 
oscillators is observed in many complex biological, chem¬ 
ical, physical and sociological systems with different ori¬ 
gins of rhythmical activity and different mechanisms of 
coupling. Depending on the context the term ’’oscilla¬ 
tor” may refer to neurons [lj, cardiac pacemaker cells 2|, 
yeast cells :3f], flashing fireflies Q, chirping crickets 5jj, 
applauding spectators 6f, nano-mechanical resonators 
0) lasers [§], superconducting Josephson junctions 0, 
etc. The concept of spontaneous synchronization has 
proved to be useful in designing of artificial networked 
systems capable of self-organization in the absence of 
any centralized control mechanism M- A common ex¬ 
ample of these systems is the so-called wireless sensor 
networks that typically consist of a number of spatially 
distributed devices with computational, sensing, memory 
and wireless communication capabilities 0- Many in¬ 
dustrial and civil applications of sensor networks require 
global clock synchronization of sensors [l2] or their dis¬ 
tributed consensus on certain quantities of interest fl3i] . 
A good synchronization scheme should provide the ro¬ 
bustness with respect to failure of an individual sensors 
and take into account the diversity of sensor’s properties 
and the presence of noise. The nature-inspired approach 
to this problem is completely distributed synchronization 
strategies based on the mutual coupling between sensors 
[14| - |l8j . In self-organizing sensor network the nodes are 
coupled through the exchange of information: each node 
transmit its own state and receive the signals from the 
neighbors. The knowledge of this information allows to 
single node to update itself occasionally in accordance 
with some correction rule so as to facilitate the global 
synchrony. Another important technological application 
of synchronization phenomenon lies in the field of smart 
power grids [IjJ |2Cj|. The modern electrical grids consist 
of thousands of generators and substations linked across 


large distances. To ensure stable operation and efficiency 
of the entire system, all the components that generate 
electricity must operate at the same frequency. In a de¬ 
centralized smart power grid, the nodes could communi¬ 
cate among each other and correct their parameters in 
order to achieve and maintain the synchronized state. 

Motivated by the discussed applications, in this work, 
we investigate an ensemble of phase oscillators for which 
the mutual coupling has a form of phase correction. Our 
model implies that each oscillator evolves autonomously 
except for discrete times when the instantaneous correc¬ 
tion events occur. We focus our attention on the simple 
correction rule which prescribe to oscillator to reset its 
phase to the ensemble-averaged value at a typical rate 
which is determined by the product of the coupling con¬ 
stant and the degree of global phase coherence. This 
basic model is inspired by the famous Kuramoto model 
in which the effective coupling force is proportional to the 
amplitude of the mean field 21]. The central result of the 
work is a synchronization transition upon the change of 
coupling constant. We demonstrate that for sufficiently 
strong coupling a partially synchronized state continu¬ 
ously bifurcates from incoherence. It is also shown how 
the synchronization properties of oscillators change in the 
case of non-linear relation between the correction rate 
and the degree of synchrony. 


II. MODEL FORMULATION 

Consider a diverse ensemble of N noisy oscillators 
which are solely characterized by the phase variables 
ipi £ [—7r,7r], where 1 < i < N. We will use the Ku- 
ramoto’s order parameter [ 2 J, 
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to describe the collective rhythm produced by the whole 
system. The magnitude 0 < r < 1 of this complex pa¬ 
rameter measures the degree of macroscopic coherence, 
while ip is the average phase. In the absence of coupling 
the dynamics of the ith oscillator is given by 

dt^Pi = O, + £,i, (2) 

where the intrinsic frequency f 2 ,; is randomly chosen from 
some probability density g(O) and £j(t) is the Gaussian 
white noise, 

<&(*)> =0, (Zi(ti)Zi(t2)) = 2D6 ij 6(t 1 -t 2 ). (3) 

The the angular brackets denote averages over statistics 
of the noise and the non-negative parameter D measures 
its intensity. Dropping the index, we introduce the one- 
oscillator probability density of a phase distribution as 
n(ip, 0, t) = (5(tp — <^(£))), where <p(t) is a solution of the 
stochastic equation © for a fixed realization of noise. 
The time evolution of this function obeys the following 
Fokker-Planck equation [23j|] 

d t n = Dd^n — Qd^n. (4) 

Note that n is nonnegative, 27r-periodic in ip, and satisfies 
the normalization condition 

+7T 

J n(ip, fl, t)dip = 1. (5) 

— TV 

In the thermodynamic limit N —> oo, the order parame¬ 
ter © can be expressed as 

+ 00 TV 

re 1 ^’= J dflg( fi) J e lip n(ip,fl,t)d(p. (6) 
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For an arbitrary initial distribution n(tp, O, 0) the so¬ 
lution of Eq. © is given by the convolution n(ip, O, t) = 
J — P>' , t)n(ip', O, 0 )d<p' , where the Green function 

1 1 00 2 

GUp,Vl,t) = - 1 — V e~ m Dt cos m(tp — fit), (7) 

27T 7r 

771=1 

obeys the Fokker-Plank equation © together with the 
initial condition G(<p, 0,0) = 5(<p). It follows from © 
that any initial state relaxes exponentially fast to the 
uniform distribution n(vs, O) = 1 / 2 - 7 T, which corresponds 
to zero order parameter ©. Predictably, the diverse 
population of non-interacting noisy units behaves inco¬ 
herently. 

Next, let us assume that oscillators are coupled and 
update itself in order to improve the degree of global syn¬ 
chrony. Namely, the ?’th oscillator evolves accordingly to 
the equation ©, but from time to time it instantaneously 
shifts the phase to the current ensemble-averaged value 
ip. We treat the updating time instants as completely 


random and statistically independent for different oscil¬ 
lators. The intensity of phase correction is characterized 
by a typical rate of updates a per oscillator which, in 
general, may change over time. Under these assump¬ 
tions, the phase dynamics can be viewed as a drifting 
diffusion process on a circle in the presence of stochastic 
resetting [24j with ip and a playing a role of resetting 
position and resetting rate, respectively. We, thus, write 
the following master equation in the limit N —> oo 

d t n = Dd^n — fld^n — an + a5(ip — ip), ( 8 ) 

where ip is determined by © and the one-oscillator prob¬ 
ability density n(<p,VL,t) implies the averaging of phase 
dynamics over both the noise and stochastic phase cor¬ 
rection. The structure of the equation © is quite trans¬ 
parent: the first and second terms in the right hand side 
correspond to a diffusion and uniform rotation at natural 
frequency (cf. with Eq. ©), while the third and fourth 
terms are related to updates and represent a negative 
probability flux out of each point ip and a corresponding 
positive probability flux into ip = ip(t). Taking into ac¬ 
count the 27r-periodicity of function n in tp it is easy to 
show the conservation of total probability ©. 

The equation © constitutes the basis of our statistical 
model for coupled oscillators. The focus of this study is 
the synchronization properties of the system in the case 
when the correction rate a is a function of the mean-held 
amplitude r. It is useful to start with the description of 
the steady-state dynamics arising from the Eq. © with 
constant parameter a. 


III. STEADY-STATE ANALYSIS 

Suppose that the rate of updates a is time- 
independent. Since the phase correction promotes the 
consistent behaviour of oscillators, the incoherent state 
n(ip,Q) = 1/(27r) does not solve the Eq. © at a ^ 0. 
If frequency distribution g(fl) has a single maximum at 
frequency flo and is symmetric, it is reasonable to seek 
n(p, 12, t) in the form of some stationary profile uniformly 
rotating at frequency Do- Then passing into the rotating 
frame ip —> p + Oo t and neglecting the time fluctuations 
of the order parameter at N —> oo one can set ip = 0 
without loss of generality. In result, the problem reduces 
to the stationary equation 

Dd^n — uidtpU — an + a5(ip) = 0, 

supplemented by the periodic condition n(— 7 r, w) = 
n(n,uj), where u> = O — Oo- It is straightforward to find 
the following solution 

! A\ exp[ 7 !<p] + A 2 exp[ 72 vs], 0 < ip < n, 

Bi exp[ 717 )] + B 2 exp[ 72 vs], -7t < ip < 0, 

(9) 
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in which 71,2 = (w ± \/w 2 + 4aD)/2D and 


A 10 — 

a 

( 10 ) 


(6271-7 !,2 - l)Vw 2 +4 OLD 

Bl.2 = 

ae 2 ^ 1 - 2 

( 11 ) 
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It is noteworthy, that for a / 0 this stationary distri¬ 
bution is nonequilibrium since phase correction produces 
ongoing circulation of probability. 

We can now calculate the order parameter © by using 
mm- Since i[> = 0 by assumption, one readily obtains 


+ OO 

r = a(a + D) j dui 


g(fl o + uj) 
lo 2 + (D + a) 2 


( 12 ) 


In a particular case of a Lorentzian distribution g(fl) = 
( 7 / 7 r )( 7 2 + (£1 — Hq) 2 ) we hnd exactly 


a 

a + D + 7 ’ 


(13) 


where P(t) is the probability density of random variable 
r. After the substitution the expressions (Ha, mi and 
if) = 0 into © we derive the following simple formula 


+ OO +OO 

r = J J dojd,Tg(uj + tto)P(T)e~ DT cosojt, ( 17 ) 

—00 0 

which allows us to readily determine the stationary co¬ 
herence r once P(t) is known. For stochastic phase cor¬ 
rection with Poisson rate of updates a we have P(t) = 
ae~ aT . Then integration of d3 over t yields exactly 

da. 

As an example, let us apply the Eq. da to the case of 
independent periodical correction processes. We assume 
that the sequence of updating events {ti(fc)}fc)Lo °f *th 
oscillator is given by U(k) = Ti + kT , where T is updating 
period and r,; is a random variable uniformly distributed 
over the interval [0,T]. Then P(t) = T~ l for 0 < r < T 
and P{t) = 0 otherwise. Performing integration in (fTTl) 
with g(fl) = (7/7t)( 7 2 + (f2 — Qq) 2 ) one obtains 


whereas for a normal distribution with standard devia¬ 
tion a one obtains 


1 _ e ~(D+l)T 

{D + l)T 


(18) 



where erfc denotes the complementary error function. 

Importantly, the result da is extracted from the equa¬ 
tion ® which is justified by the assumption that updates 
are the random Poisson events. There is a more general 
way to derive the steady-state coherence r which is ap¬ 
plicable to the case of arbitrary stationary statistics of 
updating time instants as long as these instants are inde¬ 
pendent for different oscillators. As previously, we sup¬ 
pose that density profile of oscillators in the phase circle 
is stationary in the the rotating frame, ip —> + t, and 
set ^i = 0. Then the oscillator dynamics may be thought 
of as a renewal process: each update resets the phase of 
oscillators to ip = 0 , whereas in the inter-update time 
intervals the evolution of phase is given by superposi¬ 
tion of diffusion and detuning drift. The Green function 
© rewriting in the rotating frame yields the probabil¬ 
ity density of phase distribution for a randomly chosen 
oscillator 

1 I 00 

G(<p,w,r) = - + - V e -™ 2 ^ cog m((p-uT), (15) 

Z7T 7r ' 

m= 1 

where r denotes the time elapsed since the last update 
of this oscillator. To obtain the steady-state distribution 
n(ip,oj) one should average the Eq. (fl5ll over different 
realization of the phase correction process, i.e. 

+ OO 

n(tp,u>)= J P(T)G(p,uj,T)dT , (16) 

0 


In agreement with intuition, Eqs. m3D, o and m 
give r —> 0 (incoherence) and r —> 1 (perfect synchrony) 
at a,T _1 —► 0 and a,T~ l 00 respectively. 


IV. STATE-DEPENDENT CORRECTION RATE 

We now turn to the more non-trivial case when the 
correction rate a is not fixed but depends itself on the 
state of population. Let us assume that updates occur 
the more frequently, the more higher the current degree 
of phase synchrony r. Apparently, this coupling scenario 
sets up a positive feedback so that as the population 
becomes more coherent, r grows and so the correction 
rate a increases, which tends to synchronize the oscil¬ 
lators even stronger. The situation is quite similar to 
classical Kuramoto model where the mean-field coupling 
becomes stronger with the growth of global synchrony. 
It is well-known that Kuramoto-like coupled oscillators 
exhibit spontaneous synchronization provided the inter¬ 
action constant is large enough. The major question we 
wish to address here is whether there is a similar phe¬ 
nomenon in the discussed case of coupling through the 
coherence-induced phase correction. The answer is pos¬ 
itive, below we demonstrate that for sufficiently strong 
coupling the phase correction compensates the destruc¬ 
tive effects of noises and frequencys diversity and a par¬ 
tially synchronized state emerges. We focus our attention 
on the linear relation between coherence and correction 
rate, i.e. a = kr , where k > 0 is the interaction con¬ 
stant. The nonlinear generalizations of the model will be 
discussed briefly in Section [V] 
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A. Self-consistency analysis and numerics 


One may expect that starting from an arbitrary initial 
condition the population of oscillators will settle into a 
steady state with constant coherence r and average phase 
%/} running at frequency fio- Then the correction rate a 
becomes time-independent and the results of the previous 
section can be straightforwardly applied. Replacing a by 
kr in Eq. 113, we find the following self-consistency 
condition for stationary coherence r 


+ OO 

r = kr[kr + D) j dco 


g(flp + ui) 

LO 2 + (D + kr) 2 ' 


(19) 


A trivial zero solution r = 0, corresponding to a incoher¬ 
ent state, holds for any value of k. A second brunch of 
solutions bifurcates continuously at critical coupling 



k c 


1 

D 



+ ui) 
uj 2 + D 2 


— OO 


-1 


( 20 ) 


obtained by letting r —► +0, and corresponds to a par¬ 
tially synchronized state with non-zero r. 

Remarkably, formula (TSUI) has the same structure as 
that for the critical coupling in the noisy Kuramoto 
model pl|. At the same time, the models differ by the 
critical behaviour of the order parameter. An expan¬ 
sion of the right-hand side of Eq. 03 in powers of r 
yields the linear scaling law r oc k — k c in the vicinity 
of transition. In contrast, in Kuramoto model the order 
parameter of the bifurcating branch obeys the square- 
root law r oc y/k — k c near threshold. Note also that 
how r grows with k strictly depends on the statistics 
of updates. For the sake of illustration, suppose that 
g(kl) is a Lorentzian distribution. Then, for stochastic 
phase correction we find from (1131) exactly r = {k — k c )/k 
at k > k c = D + 7. However, for periodic phase cor¬ 
rection with coherence-dependent frequency of updates 
1 /T(r) = kr, one obtains from (fl8l) the non-trivial solu¬ 
tion r = —(k/k c ) ln~ 1 (k/k c — 1) at k > k c = D + 7. 

The emergence of spontaneous synchronization is con¬ 
firmed by simulations of large but finite population of os¬ 
cillators. Our numerical scheme use temporal discretiza¬ 
tion of Eq. m with time step At so that is the ap¬ 
proximation of ipi(nAt), where n = 0 , 1 , 2 ,.... Then, the 
complex order parameter is r n e z ^ n = (1/AO )Cj=i ■ 
The random forces are modeled by the telegraph pro¬ 
cesses whose correlation times are equal to time step du¬ 
ration At. Specifically, the values of 0: inside the nth time 
step are chosen to be independent random constants 
with identical normal distributions. To ensure a given 
value of the diffusion coefficient D , one should choose 
(£r 2 ) = 211/At. At each time step we generate the set 
of N random numbers r", r|*,..., t]J- having distribution 
a n e~ anT with a n = kr n . If r” < At, we switch the phase 
of tth oscillator to the current ensemble-averaged value 


FIG. 1. The degree of phase coherence r in dependence on 
coupling strength k for population of N noisy oscillators with 
normal distribution g(Sl) of natural frequencies. The solid line 
shows the solution of the self-consistency equation Cl) corre¬ 
sponding to thermodynamic limit N —> 00. Numerical results 
are shown as points, and are taken from the simulations of 
N = 1000 oscillators. The simulations were performed for 
2 x 10 5 time steps, starting from uniform initial distribution 
of phases over the interval [—7r, n]. The values presented in 
the figure are averages over the last 1.8 x 10 5 steps. 


ij) n . Thus, the rules for the evolution of the system are 
as follows 

f < + (Vi + C)A t, if > At, 

^ +1 = { 

{ V’n + O^ + erXAt-T"), if rf < At. 

In numerics the initial phases of N = 1000 oscillators 
were uniformly distributed over the interval [— n, 7r] and 
the discrete time step At = 0.001 was used. The diffusiv- 
ity was D = 1 and the natural frequencies fli were cho¬ 
sen from normal distribution g(ft) with expected value 
= 0 and standard deviation a = 0.5. Figure [T| repre¬ 
sents numerical results for the time-averaged coherence r 
as a function of coupling constant k in comparison with 
theoretical curve for infinite-A system. It can be seen 
that level of synchrony in population remains low until 
k reaches a critical value k c ~ 1.2, above which r mono- 
tonically increases towards its asymptotic value of unity. 
The non-zero values of r below k c in the simulation reflect 
fluctuations due to finite N. 


B. Global dynamics in homogeneous system 

The self-consistent arguments presented above do not 
provide any information about the stability of incoherent 
and partially synchronized states. The direct analysis 
of the stability properties is possible when the diversity 
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of natural frequencies is negligible, i.e. g(Ll) = 6(fl — 
f2 0 ). Then, using Eq. © one can derive the following 
closed-form equations associated with angular and radial 
motions of the order parameter re*^ = e ltp n{<p,t)dp 


= fto, 

r = —Dr — ar + a. 


( 21 ) 

( 22 ) 


The equation (l2ll) suggests that the average phase just 
uniformly rotates at frequency fio- Substituting a = kr 
we find easily the solution of amplitude equation (l22l) 


r(t) 


(k - D)e ( - k ~ D '> t r 0 
k — D — k(l — e < ' k ~ D '> t )r 0 ’ 


(23) 


where r 0 = r(0) is an initial condition. For k < D the 
coherence r always decays to zero as t —> oo, but for 
k > D incoherent state loses stability and one obtains 
another attracting point r = (k — D)/k corresponding to 
partially synchronized state. Since g(Ll) = <5(S2 — flo) by 
assumption, the critical value k c = D is consistent with 

If the natural frequencies of oscillators are non¬ 
identical, there is no closed description of global dynam¬ 
ics in terms of the order parameter. In the next subsec¬ 
tion we propose perturbative approach which allows us 
to demonstrate that incoherent solution becomes linearly 
unstable as the strength of the coupling goes beyond the 
threshold ([20]) . 


C. Linear stability analysis of incoherence 

In what follows we look at the oscillator dynamics in 
the frame moving with frequency Ho. Let us consider a 
small perturbation of the uniform incoherent state 

n(<p,u,t) = ^ +ep(ip,u,t), (24) 

where eCl. At lowest (linear) order in e the evolution 
of perturbation is governed by equation 

dtp = Dd^p- ud v p- + kr5(p - ip). (25) 

To proceed we write a 27r-periodic function p as a super¬ 
position of Fourier harmonics 

OO OO 

p{<p,u,t) = £ c m (w, t)e imv + £ (26) 

m =1 m =1 

Note that normalization condition ([5]) automatically pro¬ 
vides Co = Cq = 0. When (l26l) is substituted into (|25|) , 
we obtain the following equation for the amplitude c m 

dtc m = — ( m 2 D + imui)c n + (27) 

The evolution equation for cj^ is just the conjugate of 

(EH). 


Next, substituting (l26l) into (HD yields 


+oo 




= 27T 


diog(io + flo)ci(w, t). 


(28) 


Thus, only the first harmonic of the Fourier decompo¬ 
sition, which is the so-called fundamental mode, con¬ 
tributes to the order parameter. 

From EH) and C51) one easily sees, that equation for 
the amplitude c\ (and c^) is uncoupled 

+oo 

<9tCi = — (D + iuj)ci + k J dvg(v + flo)ci(^, t). (29) 

— OO 


Let ci(ui,t) = a(w)e At , then we obtain La = A a, where 
L = — (D + ioj) + k dvg(v + fl 0 )- The spectrum of 

the operator L was constructed in [26} in the context of 
noisy Kuramoto model jH}. They show that the con¬ 
tinuous part of spectrum lies in left half-plane for any k , 
thus corresponding modes never cause instability. In con¬ 
trast, the location of discrete spectrum depends strongly 
on interaction constant k. When k exceeds a threshold 
m the discrete real eigenvalue A > 0 appears so that 
incoherent state becomes unstable. 

Now we turn to the analysis of high order amplitudes. 
The equation (l27l) may be solved easily in terms of r(t), 
ip(t) and the initial condition c m (w, 0). The result is 

CmM = c m (w, 0)e-< m2D+im “>* + 

t 

-L-h e -( rn2D + irnuJ ) t 

0 




The order parameter is completely determined by the 
fundamental mode through Eq. ([28]) . thus we put 
r(t)e*'^( t ^ = r(0)e i ^(°) + ' u , where the eigenvalue A belongs 
to the spectrum described above. After integration one 
obtains 


Cm (w 1 1) 

T ( o m (cc, 0) 


kr( 0) 


-,— 2771 - 0 ( 0 ) 




m 2 D + A' + im(u> — A") 
kr{ Q) e ~ im ^) 


m 2 D + A' + im(uj — A") 


— (m 2 D-\-i 




where A 7 and X" are the real and imaginary parts of A, 
respectively. Notice that the denominator in (1301) never 
vanishes since m > 1 and A' > — D for any k (26[. The 
second term in the rhs of rapidly dies out, while 
the behaviour of the first one depends on the coupling 
strength k. At k <k c the spectrum A lies in the left half¬ 
plane so that c m (t) always decays to zero. For k > kk 
the aforementioned positive discrete eigenvalue provides 
the exponential growth of 

Let us summarize the results of this subsection. The 
incoherent state goes linearly unstable to perturbation 
involving first harmonic for k > k c , where the critical 
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FIG. 2. The phase coherence r in dependence on the coupling 
strength k for infinite ensemble of identical noisy oscillators 
with nonlinear phase correction rate a = kr” 1 , where m — 2 
(left subplot) and 1/2 (right subplot). The diffusion coeffi¬ 
cient is D = 1. The solid and the dashed lines correspond to 
the stable and unstable branches respectively. 

coupling k c coincides with the prediction m °f se if~ 
consistency analysis. Above the threshold the instability 
of the first harmonic induces also the growth of the high 
order amplitudes. 


V. NONLINEAR GENERALIZATION 

So far, we have considered the linear relation between 
correction rate a and amplitude of the order parameter r. 
One may expect more complex and interesting collective 
behaviour in system with non-linear dependence a(r). 
Let us discuss as an illustration the quadratic model 
a = kr 2 focusing on the case of identical oscillators, i.e. 
g{ 0) = i5(fl — Slo). Then, the amplitude equation E21) 
reads: 

f = — Dr — kr 3 + kr 2 . (31) 

At k < 4 D the only time-independent solution is r = 0 
and r(t) decays to zero as t —> oo for any initial con¬ 
dition r(0) = ro- However, for k > 4D there are three 
fixed points: 0 and r± = (1 ± sj\ — 4D/k)/2. Thus, 
two partially synchronized branches bifurcate discontin- 
uously from r = 1/2 at k c = 4 D, as it is shown in Fig. 
[2j The incoherence (r = 0) and the upper synchronous 
state (r = r+) are stable, and the lower synchronous one 
(r = r_) is unstable. This means that the initial con¬ 
figurations with ro < r_ relax to incoherence as t —f oo, 
while those with ro > r_ converge to the synchronized 
state at long-times. The described bistability is clearly 
observed in numerical simulations of finite-IV system, see 
Fig. [3l Note also that in more general case a = kr m with 


FIG. 3. Temporal evolution of coherence r(f) in population of 
N = 1000 identical oscillators with quadratic phase correction 
rate, a = kr 2 . The diffusion coefficient is D = 1 and the 
coupling strength is k = 5 > k c . Different curves stand for 
different initial configurations of phases. The upper and the 
lower dashed lines correspond to stable and unstable partially 
synchronized states of infinite-IV population respectively, see 
Sect. \V\ 

m > 1 the critical coupling is k c = Dm m /(m — l) m_1 
and two partially synchronized branches bifurcate from 
r = ( in — l)/m. Here again the system is bistable with 
respect to initial conditions. 

For the second particular example a = kr 1//2 the am¬ 
plitude equation (E^l) yields 

r = —Dr — kr 3//2 + kr 1 ^ 2 . (32) 

The incoherent solution r = 0 is unstable, while the only 
non-trivial solution r = [( \/D 2 + 4k 2 —D)/(2k)] 2 is stable 
and exists for any k. Thus infinitely small coupling leads 
to spontaneous synchronization of system. This feature 
is common for all models a = kr m with 0 < m < 1. 

As it was stated in the previous section, if oscillators 
have different natural frequencies, there is no closed evo¬ 
lution equation for the order parameter. Nevertheless, 
the stationary values of r can be found by substitution 
a(r) into (USD and solving the resulting self-consistency 
equation. It is easy to see from m that when the fre¬ 
quency diversity is described by a Lorentzian distribu¬ 
tion g(Q) one should just replace D by D + 7 in the 
above formulas for identical oscillators. Notice, however, 
that these arguments do not indicate which of the found 
branches are stable. 


VI. CONCLUSION 

In this work, we proposed and investigated a new 
model of coupled non-identical noisy phase oscillators. 
We assumed that the mutual coupling manifests itself as 
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a tendency of each oscillator to adjust its phase to the 
ensemble-averaged value through instantaneous phase 
shifts which occur at random updating time instants. 
Then in the thermodynamic limit the time evolution of 
the one-oscillator probability density function is governed 
by the equation (j4j. The key ingredient of our model is 
the coherence-dependent rate of updates: the more pro¬ 
nounced the coherence of the current state of population, 
the more frequently the phase shifts occur. We mainly 
concentrated on the simplest linear case when the typ¬ 
ical rate of updates is given by the coupling constant 
multiplied by the coherence. The distribution of natural 
frequencies of oscillators was assumed to be unimodal 
and symmetric. Our study revealed that the system 
exhibits the transition to spontaneous synchronization 
in a fashion similar to the mean-field Kuramoto model. 
Specifically, using the self-consistency arguments for an 
infinitely large population, we evaluated analytically the 
critical coupling strength (EU1) that marks the continu¬ 
ous onset of a partially synchronized solution. The same 
threshold follows from the linear stability analysis of in¬ 
coherent state. In the vicinity of the critical value, the 
linear scaling of the order parameter is found. For the 
sake of illustration, we performed the numerical simula¬ 
tions of large number of oscillators undergoing coherence- 
induced phase correction. The analytical predictions for 
the degree of phase coherence in dependence on the cou¬ 
pling strength turned out to be in good agreement with 
numerical results, see Fig. ©■ 

For homogeneous population of oscillators, we found 
a closed description of global dynamics, which consists 
of two uncoupled evolution equations ET1) and (l22l) for 
the amplitude and phase of the complex order parameter. 
Based on this result, we were able to demonstrate easily 
the global stability of a partially synchronized state in 
the system of identical oscillators. 


Finally we considered the direct generalization of the 
model to the case of a nonlinear coupling. When the re¬ 
lation between the correction rate and the coherence is 
given by a power law with exponent larger than unity, 
the synchronization transition is discontinuous (Fig. [2]) 
and above the threshold the oscillators exhibit bistability 
of global activity between incoherence and a partial syn¬ 
chronization (Fig. [3j). If the correction rate grows with 
the coherence slower than linearly, the system shows the 
non-zero level of synchrony for arbitrary small coupling 
constant (Fig. [2]). 

The model presented here can be further extended in 
many ways. Obviously, the knowledge of the current 
states of nearby nodes cannot be assumed perfect in a 
networked applications because of communication and 
processing time delays and information losses. From the 
perspective of the proposed theoretical framework this 
means that the delta-function in the rhs of our basic 
equation © should be replaced by some spread func¬ 
tion which depends on the previous states of the sys¬ 
tem. Another natural extension of the model includes 
the short-range interaction effects for regular and com¬ 
plex networks. In future studies, it is important also 
to consider the synchronization properties of system in 
which the updating time instants of different oscillators 
are correlated. 

Here we have restricted ourselves to the case of ’’at¬ 
tractive” phase correction which attempts to mutually 
synchronize the oscillators. There are many practical sit¬ 
uations in which the formation of a coherent state in en¬ 
semble of oscillatory units is unacceptable and should be 
avoided [13,HI]. One may expect that ’’repulsive” phase 
correction (i.e. the phase shifts to ip+ir rather than to ip) 
is a good strategy to break the undesired synchronization 
of Kuramoto-like coupled oscillators . 
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